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We study pairing vibrations in 1^,20,22 q 42,44,46 nuclei solving the time-dependent Hartree- 
Fock-Bogoliubov equation in coordinate space with spherical symmetry. We use the SLy4 Skyrme 
functional in the normal part of the energy density functional and a local density dependent func- 
tional in its pairing part. Pairing vibrations are excited by two-neutron transfer operators. Strength 
distributions are obtained using the Fourier transform of the time-dependent response of two-neutron 
pair-transfer observables in the linear regime. Results are in overall agreement with quasiparticle 
random phase approximation calculations for Oxygen isotopes, though differences appear when in- 
creasing the neutron number. Both low lying pairing modes and giant pairing vibrations (GPV) are 
discussed. The GPV is observed in the Oxygen but not in the Calcium isotopes. 

PACS numbers: 21.60.Jz, 24.30.Cz, 21.10.Re, 25.60.Je 



I. INTRODUCTION 

Energy density functional (EDF) approaches like the 
Skyrme-Hartree-Fock model [l], have proved to be suc- 
cessful to describe bulk properties of nuclei over the nu- 
clear chart 0]. Recent computer improvements allow 
large scale EDF calculations of nuclear structure [3, Q 
and reactions [1, 0]- In these approaches, one restricts 
the many-body wave- functions to a subset of the Hilbert 
space on the one hand and guess a nuclear EDF on the 
other hand. A commonly used technique is to break sym- 
metries to enrich the variational sub-space and improve 
the description of nuclear structure. As an example, 
breaking gauge invariance associated to particle number 
conservation yields the Hartree-Fock-Bogoliubov (HFB) 
formalism 0, [^, [l^. This technique allows for the de- 
scription of superfluidity in ground states of open-shell 
nuclei. 

Extensions to treat collective excitations in the pres- 
ence of pairing correlations are possible in the frame- 
work of the Quasiparticle Random Phase Approximation 
(QRPA) which has been widely used in nuclear struc- 
ture studies dH, m [H, H, m, m m. This approach 
and its zero pairing counterpart (RPA) give reasonable 
estimates of giant resonances though improvements are 
necessary to reproduce fine structures [l8|. In fact, the 
QRPA can be obtained from the linearization of the time- 
dependent Hartree-Fock-Bogoliubov (TDHFB) equation 
which provides a self-consistent evolution of an indepen- 
dent quasiparticles state. There is a consistency require- 
ment that QRPA and the static limit of TDHFB should 
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use the same effective interaction. This is a natural fea- 
ture of TDHFB that the same EDF can be easily used 
in the static and dynamical calculations thanks to the 
structure of the TDHFB equation. This is not always 
the case in (Q)RPA calculations where spin-orbit and 
Coulomb parts of the residual interaction are often omit- 
ted (see discussion in [l^ and references therein), which 
may affect collective modes !15, 20, 21, 22]. 

Unlike in condensed matter where, for instance, TD- 
HFB has been applied to study dynamics of Bose- 
Einstein condensates [H, [13] , explicit time evolutions of 
nuclei including pairing are sparse and usually limited to 
the BCS ansatz of superfluidity with simple functionals 
[25I [2^ [27} or to simple systems 2&1 . Only recently, a 
numerical method of solving TDHFB with the Gogny in- 
teraction [2^ has been proposed to study quadrupole 
oscillations using a harmonic oscillator basis [3(j ]. 

At the limit where pairing is neglected, however, ex- 
tensive calculations of nuclear dynamics have been per- 
formed using the time-dependent Hartree-Fock (TDHF) 
formalism introduced by Dirac in 1930 [sij . In this ap- 
proach, one considers the dynamics of independent par- 
ticles in a self-consistent mean-field generated by all the 
others. The use of Skyrme EDF 32] allowed recent realis- 
tic calculations of both collision mechanisms d, [tI. [33l . Is^ . 
m, HE [13, [11] and giant resonances [H, [H, [Mll^ For 
instance, TDHF has been used to study anharmonicities 
in collective motion which are beyond the RPA range of 
applications [H, \^ . It is then an appealing challenge to 
repeat these works using TDHFB to investigate the role 
of pairing correlations in nuclear dynamics. In particular, 
"How do they affect low-energy reaction mechanisms ?" 
is still an open question. However, full three-dimensional 
TDHFB codes for collisions are still prohibitive at mo- 
ment. For nuclear structure purposes, pairing vibrations 
and rotations [43, WJ, W3 also demand theoretical inves- 
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tigations to reach realistic predictions [4g| . Of particular 
interest are high-lying pairing collective modes like gi- 
ant pairing vibrations (GPV) [13]. These modes can be 
viewed as coherent sums of two-quasiparticle excitations 
across a major shell. They are probed by two-nucleon 
transfer reactions [4l,|4§|. The GPVs arc still unobserved 
experimentally but they are experiencing a renewed inter- 
est since recent theoretical developments predicted that 
the use of radioactive ion beams could provide better 
conditions for their studies [13, HH US] ■ 

In this article, we solve for the first time the TDHFB 
equation with a full Skyrme functional in the normal part 
of the EDF and a local density dependent one in its pair- 
ing part. As a first application, we investigate pairing 
vibrations using the linear response theory. Our model 
is then equivalent to a fully consistent QRPA including 
spin-orbit and Coulomb interactions also in the particle- 
hole channel. In scction|TTl we recall the TDHFB formal- 
ism and the choice of the EDF. We present in section Hill 
numerical implementations in spherical symmetry. Then, 
we discuss in section IIVI our choice of observables asso- 
ciated to pairing vibrations in the framework of the lin- 
ear response theory. Finally, we present the results for 
18,20,220 and 42.44,46ca isotopcs in section |V] before to 
conclude in section IVll 



II. FORMALISM 



A. TDHFB equation 



The TDHFB equation can be derived starting from the 
action between an initial and final time U and tf 

s = J\t{^mh^^-H\^{t)) (1) 

and writing the variational principle 55 = in the sub- 
space of quasiparticle vacua. For each state \'^) of this 
sub-space, one can find a basis of quasiparticle annihila- 
tors {/?} such that $^\'^) = for all /i [M]- The latter can 
be related to the particle creation and annihilation oper- 
ators {a^ , a} through the Bogoliubov transformation fTo| 



E 



(2) 



where the matrices U and V are such that the quasi- 
particle operators {/3^,/3} fulfill the canonical anti- 
commutation rules for fermions. 

The variational principle leads to the TDHFB equa- 
tion [si 



ot 



The generalized one-body density matrix reads 

7^(^) = 



p{t) K{t) 

-n*{t) l-p*{t) 



(3) 



(4) 



where p^^, = (Vl^t'^Ml''/') matrix elements of the 

normal density and = {ij\ava^\4') are the elements of 
the pairing tensor. The generalized one-body Hamilto- 
nian H, has a block structure and can be writen in terms 
of the Hartree-Fock (HF) Hamiltonian h and the pairing 
field A 



h 

-A* 



A 

-h* 



where 



S£[p, K, K* 
Spi^f, 



and A„^ = 



S£[p, K, K*] 



6k 



(5) 



(6) 
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In the above HFB formalism, the functional £[p, k, k*] is 
the expectation value of the exact Hamiltonian H on the 
quasiparticle vacuum 

A more practical form of the TDHFB equation is found 
by recasting these equations in terms of the quasiparticle 
components U and V introduced in Eq. ([2]) 



at V 



E 



-A* -h 



(7) 



B. EDF approach 

In nuclear physics, the above formalism should be mod- 
ified, in particular to take into account the short range 
repulsive part of the nuclear interaction which makes 
the mean-field HFB approach irrelevant with the bare 
interaction. This is a reason why one generally re- 
places £[p,K,K*] by an effective EDF fitted on nuclear 
properties without invoking directly the underlying ex- 
act Hamiltonian. Moreover, in the spirit of the density 
functional theory ^5^, [s^ , this procedure allows to 
include many-body correlations. 

Let us decompose the total energy into kinetic, 
Skyrme, Coulomb and pairing parts (see appendix [XI for 
an explicit expression of each component) 



kin 



£ 



Sk 



£ 



Coul 



£. 



■pair • 
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We choose the SLy4 parameterization [57] of the Skyrme 
functional 's^l including time-odd densities [58|. The 
Coulomb energy includes direct and exchange terms. The 
latter is estimated using the Slater approximation [59| . 
The three first terms of Eq. ([5]) depend only on the nor- 
mal densities. It is convenient to express the pairing en- 
ergy £pair using the anomalous density 

Pq{rs, r's') ^ -2s'Kg(rs, r' - s') (9) 

where s the projection of the spin and q the isospin. We 
use a local pairing functional (see, e.g., |2] and references 
therein) 



-'pair 
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where pq(r) = J2s Pqi^^^^^^) and p,(r) = EsP9(rs,rs) 
are the local parts of the normal and anomalous densi- 
ties with isospin q respectively, po(r) = J^qPqi'"^) is the 
scalar- isoscalar density and g is the pairing coupling con- 
stant. The parameters pc and 7 are adjusted to generate 
pairing correlations preferably at the surface and/or in 
the bulk of the nucleus. Such a pairing scheme yields 
pairs of nucleons of the same isospin coupled to angular 
momentum zero. The simplicity of such an EDF made 
systematic three-dimensional HFB calculations over the 
whole nuclear chart possible 0, 0] . 

However, one has to face the divergence of local pair- 
ing densities [60| . It is then necessary either to regularize 
the equations by introducing a cutoff in the quasiparticle 
spectrum or eventually to perform a more complex renor- 
malization scheme (for an example based on the Thomas- 
Fermi approximation, see [6ll. [63|). 

In our calculations, we use a cutoff to regularize the 
TDHFB equation in a quasiparticle energy window of 
80 MeV. This value allows two-quasiparticle excitations 
up to 160 MeV. The pairing parameters are fitted to re- 
produce a neutron spectral gap of 1.25 MeV in ^^"Sn, the 
pairing acting both in the bulk and at the surface of the 
nucleus (63,] . The obtained pairing coupling constant is 
g = -275.25 MeV with the parameters pc = 0.32 fm"^ 
and 7=1. 



C. Particle number conservation 

Due to the broken U{1) gauge invariance associated 
to the particle number conservation, the HFB states are 
not eigenstate of the particle number operator N. In 
static calculations, one adds a Lagrange multiplier A, 
interpreted as a chemical potential, in order to fix the 
number of particles in average. In TDHFB dynamical 
simulations, however, the particle number obeys to 



ih^{N) ^ Tr ( kA* -Ak*). 



(11) 



The definitions of A and p in Eqs. ^ and ^ respec- 
tively, together with the choice of the pairing functional 
in Eq. pO|) ensure that the right hand side of Eq. (fTTI) 
vanishes. As a consequence, we do not need to enforce 
the conservation of the average number of particles with 
a chemical potential in the TDHFB equation. 

Nevertheless, dropping this constraint in the dynamics 
induces a rotation of the Bogoliubov vacuum in gauge 
space [g^. In particular, the anomalous density of a sta- 
tionary state will carry a phase exp {—2iXt/h). Then, the 
ground state expectation values of observables which are 
linear in the anomalous density will evolve in time. For 
instance, this is the case of the observable we use to study 
pairing vibrations (see Eq. (|^^ ). The time- Fourier anal- 
ysis of such an observable in the linear response theory 
(see section irv|) will then contain a spurious peak at an 
energy Tiuj = 2X. This is a manifestation of the Goldstone 
mode due to the broken symmetry associated to particle 



number conservation. In QRPA calculations, it induces 
a spurious mode at zero energy. In order to avoid such a 
spurious mode in the linear response of TDHFB, we have 
to keep the static ground state chemical potential in the 
particle-hole field during the evolution. Finally, we note 
that the chemical potential A is easy to compute only 
for nuclei with pairing. Thus we do not consider doubly 
magic nuclei as initial states in the present work. 



III. NUMERICAL IMPLEMENTATION OF 
TDHFB 

A. Spherical symmetry 

The pairing functional defined in Eq. (jlOp couples only 
nucleons of the same isospin. We can then focus on 
semi-magical nuclei for which the spherical assumption 
is a good approximation. As a consequence, we solve the 
TDHFB equation using spherical symmetry. 

Let us recast the problem with purely local fields in 
space, spin and isospin using this symmetry. The total 
many-body wave- function being rotational invariant, it is 
convenient to write the Bogoliubov transformation in the 
spherical basis using the standard notation for quantum 
numbers n, I, j and m (we omit the isospin q in the 
notation for simplicity) 



klj7n> 



.(12) 



This definition ensures that the component m of these 
quasiparticle operators transforms under rotation as a 
tensor of rank j (see Eq. (8.79) of Ref. [H). 

We choose to solve the TDHFB equation in coordinate 
space using quasiparticle wave functions Uv{ys) = Urs,v 
and Vi,{ys) = Vrs.v defined as components of quasiparti- 
cle spinors 



Unljra (rs) 
Vnljm (rs) 



Ijm I 
'nljm I 



(13) 



where |— ) is the particle vacuum. The standard 
decomposition of single particle orbitals in spheri- 
cal coordinates writes {rs\nljm) — Rnij{r) Clijmsif^4')- 
The angular part is expressed in terms of Clebsch- 
Gordan coefficients and spherical harmonics flijms = 
{l{m — s)^s\jm)Yi^^'^ . Defining the radial quasipar- 
ticle wave functions Unji{r) = Tj2k^kn"^^^kij{r) and 
v^.iir) = (-1)'+V^,v(f")i?^.,^.(r), and using the 
property = -2s(-l)"+'-Jr!,,„_,, Eq. ^ be- 

comes 



Unljm (rs) 



1 / Unij (r) nijrr, 



FoUowing the same way as Dobaczewski et al. for the 
static HFB problem [65j . we introduce the anomalous 
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field hq{rs,r' s') = s^f*^{rs 'f's\ "^here only the pairing en- 
ergy in Eq. (fTO|) contributes. The EDF considered here 
contains only local densities (see appendixlA]). Therefore, 
the HF and anomalous fields are also local in space. Fi- 
nally, it is possible to recast the TDHFB equation ([7]) as 
a set of Schroedinger like equations for the quasiparticle 
radial wave functions u and v 

^dt \ Vnlj h* -hl^ +>^)\ ^nlj ) ^ ' 

where A is the chemical potential (see section Hi CI) . The 
expressions of the fields h{r) and h{r) and those of the 
various densities entering the EDF solved in spherical 
symmetry can be found in appendix [K\ 

B. Computational details 

The initial condition is obtained with the hfbrad 
code [1^ which solves the static HFB equation in spher- 
ical symmetry. We have constructed a time-dependent 
version of this code to solve the TDHFB equation with 
the functional described in section Hi B I The set of equa- 
tions ([H]) is solved iteratively using a one-step predictor- 
corrector method j67| and a truncation of the time prop- 
agator 

U{t,t + 6t)^exp(^^^n{t + 6t/2)^ (15) 

at 4th order in St. 

Spatial derivatives are calculated in a discretized r- 
space using seven-points formula. The numerical accu- 
racy of this method decreases increasing quasi-particle 
energy. As a consequence, this approximate derivation 
formula induces a small periodic variation of the HFB 
ground state (of the order of few tens of keV) due to high 
energy quasiparticles. We checked that this numerical ar- 
tifact disappears if we develop the wave-functions on a 
constant step Lagrange mesh j6^]. However, the latter 
method increases the numerical effort. The amplitude 
of these variations being linked to the mesh discretiza- 
tion, a mesh step of 0.15 fm has been found to ensure 
a good numerical precision with a reasonable computa- 
tional effort. In particular, the energy is conserved up to 
15 keV and deviations of the total number of particles are 
of the order of 10~^ in the present calculations. Though 
this latter value is one order of magnitude higher than 
in the TDHF case with the same numerical conditions, 
it is small enough to leave the observables of interest un- 
affected. Moreover, to avoid any unphysical contribution 
due to the approximate spatial derivative in the evolution 
of observables, we subtract from their expectation values 
the one obtained without external field. We checked that 
this procedure does not affect the physical content of the 
spectra presented in Sec. E) 

Let us precise that we use hard box boundary condi- 
tions. The latter are not optimized for a proper treat- 
ment of the continuum because they lead to a discretized 



quasiparticle spectrum. In addition, particles which are 
refiected on the boundaries of the box may interact with 
the nucleus and induce unphysical effects on the evolution 
of observables. This problem has been tackled in TDHF 
with the help of absorbing boundary conditions ^IS , ,69| . 
However, we did not found this technique to be appropri- 
ate to treat the TDHFB continuum because of the non- 
vanishing asymptotic nature of the upper component of 
the quasiparticle wave-functions in Eq. (fT3)) [65[ . Though 
needed to refine description of unbound states [l^l, fur- 
ther improvements of the boundary conditions are be- 
yond the scope of this exploratory work. 

We consider nuclei with magic proton numbers and ex- 
citations acting on neutrons only. Then the calculations 
include pairing for neutrons only. The local pairing func- 
tional couples to high angular momenta and high energy 
quasiparticle states up to the energy cutoff of 80 MeV. 
In the calculations presented here, convergence of the 
static solutions have been obtained with a maximum to- 
tal angular momentum for neutrons, for which the Bo- 
goliubov transformation is achieved, of jc = 19/2 (23/2) 
for Oxygen (Calcium) isotopes respectively. For Oxygen 
isotopes, a box radius 22.5 fm was used, in order to be in 
conditions as close as possible to the discrete QRPA re- 
sults of Ref . (E^l • Calculations for Calcium isotopes have 
been performed in a box of radius 30 fm. 
Both the mesh step, through the derivation formulae, 
and the maximum angular momentum, through the cen- 
trifugal part of the kinetic operator, constrain very much 
the time step for which the calculations are stable. The 
adopted time step used in these calculations are 0.003 
and 0.002 fm/c for Oxygen and Calcium isotopes respec- 
tively. 



IV. LINEAR RESPONSE FRAMEWORK FOR 
PAIRING EXCITATIONS 

A. Linear response theory 

The linear response theory has been widely used with 
TDHF to st udy collective vibrations in nuclei 19, 39, 4^, 
EH, m, H IH H M, in this theory, one 

computes the time evolution of an observable 

AQ{t) = {^it)\Q\^{t)) - {0\Q\0) (16) 

after an excitation induced by a small external potential 
Vextit) — £F^{t) on the ground state |0) of the system. 
The parameter e quantifies the intensity of the excitation. 
It has to be small enough to ensure the linear regime, 
i.e., the amplitude AQmax of the response must be pro- 
portional to e. In this study, the time dependence of 
the external potential is chosen to be a Dirac function 
^(t) = S{t). The excitation is then equivalent to a boost 
applied on the ground state at the initial time 

|5-(0)) = e-"^/''|0). (17) 
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The response AQ{t) to this excitation can be decom- 
posed into various frequencies lo using 

Rq{u;) = dt AQ{t) sm{ujt). (18) 



In the particular case where the operators used for the 
excitation and the observations are the same, i.e., F — Q, 
Eq. (fT5)) gives the strength distribution 

RFito) = K^-^cr) (19) 



where (a|i^|0) is the transition ampHtude between the 
ground state and the eigenstate |a) of the Hamiltonian 
and hu)a = Ea — Eq is their energy difference. 

When the excitation generates a transition to neigh- 
boring nuclei, Eq and E^ are the energies of the ground 
and excited states of the final nucleus if one keeps the 
static chemical potential in the Hamiltonian (see section 
10.1 of Ref. [53]). For instance, in the case of an addition 
of two nucleons on a nucleus of A nucleons, the energy 



of the mode reads tiujr, 



E. 



{A+2) 




where the 



ground state energy of the final nucleus is approximated 
by = E^^'^''+2X. Note that this energy may differ 



from the one obtained in a HFB calculation in the A + 2 
nucleus because of a possible rearrangement of the HFB 
field. 



B. Application to pairing vibrations 



To preserve spherical symmetry, we focus on monopole 
pairing modes, requiring a radial dependence only, i.e., 
/(r) = /(r). We choose a Fermi-Dirac spatial distribu- 
tion /(r) = (1 -I- exp ( j ,^'" )) where the parameters 
Rc = (1.27 A^l^ + 4) fm and d ^ 0.5 fm are chosen 
to allow for pair transfer on the whole nucleus on the 
one hand, and to remove unphysical high energy modes 
associated to pair creation outside of the nucleus on the 
other hand. Finally, this excitation may change the num- 
ber of neutrons at the initial time. However, deviations 
are small in the present calculations (~ 10"^ neutrons). 

We choose to follow the excitation operator F itself 
to get its strength distribution defined in Eq. p9)) . We 
also decompose the excitation operator F = J^i Ei into 
components of single particle angular momentum I 



El = J2 Jdvfir){nljm\rl){n'lj{- 



-m)\r t) 



nn jm 

^nljm'^n'lj( — m) 



xdl,„-_aL,,, „,^ + h.c. 



(23) 



where h.c. denotes the Hermitean conjugate of the entire 
expression. Computing the response of the Fi helps us 
interpret the spectra in terms of specific quasiparticle 
excitations. 

We perform the calculations over a total time interval 
of T = 1200 fm/c. In order to minimize the effects of 
the time gate on the Fourier transforms, we follow the 
protocol given in Ref. (69j . multiplying the observables 
by a time filter cos^ (w)' This procedure induces an 
additional width of 1 MeV. 



Pairing vibrations of quantum numbers 0+ can be ex- 
cited by two-nuclcon transfer reactions and have been 
studied in the small amplitude limit within the QRPA 
framework |52||. A Hermitean pair-transfer operator is 
given by [45| 

where v denotes the time-reversed state of v. 

In this paper, we consider local excitations acting on 
neutrons only and, for the sake of simple notations, we 
do not write explicitly the isospin quantum number in 
the following. The pair-transfer operator then writes in 
coordinate space 




In this particular choice, the spatial distribution /(r) is 
real. Using Eq. ([S]), the expectation value of F simply 
writes 

mWm)) - \ Jdr /(r) ipoir;t) + p*ir;t)) (22) 
where po^Y^qh- 



V. RESULTS 

In this section, two-neutron transfer in several nuclei 
is studied. To illustrate the method described in the 
previous section, we first detail the analysis in ^^O. Then, 
we present the results on neutron-rich Oxygen isotopes 
and on /-shell Calcium isotopes. 

A. Detailed analysis on 

We apply the boost of Eq. (fT7| in the linear regime on 
the HFB ground state of ^'^O. The pair-transfer operator 
F is defined in Eq. pT|) . The variation of the expectation 
value of F, obtained from Eq. (j^^ . is plotted in Fig. [1] as 
a function of time. We observe a complex evolution due 
the excitation of several modes at different energies. We 
see in the inset that a strong variation of AF(i) occurs at 
early times because all modes arc initially in phase [39] . 

The evolution in Fig.[T]is used to compute the strength 
distribution of F according to Eq. and using the time 
filter procedure described in the previous section. We 
have controlled that the extracted strength is indepen- 
dent of the excitation amplitude e. The resulting spec- 
trum is shown in Fig. dja) in solid line. We see two 
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FIG. 1: Evolution of AF{t) after a pair transfer type exci- 
tation on ^*0. The inset shows the same quantity at early 
times. 



DC 



LL 

DC 








10 20 

ha (MeV) 



30 



FIG. 2: Decomposition of the responses (in arbitrary units) 
into frequencies ui for a two-neutron pair transfer excitation 
in ^*0. (a) Strength distribution of F obtained from TDHFB 
(solid line) and the unperturbed approximation (dotted line) . 
The arrows indicate the IP3/2 (solid) and lpi/2 (dotted) deep 
hole states, (b) TDHFB responses of Fi with / = 0, 1, 2 and 
3. 



separated peaks at 6.5 and 14.3 MeV and several peaks 
between 20 and 26 MeV. 

To understand the effect of the residual interaction 
which is included in TDHFB, we have computed the so- 
called unperturbed response to the pair transfer excita- 
tion. The latter is equivalent to Eq. (fT9|) if one assumes 
that the states \a) are two-quasiparticle excitations of 
the type Ifiv) = ,5j/3^|0) where $1 creates a quasiparticle 



eigenstate of the static HFB Hamiltonian on its ground 
state |0). In this approximation, the energy of the tran- 
sition to the state Ifiv) is the sum of the quasiparticle 
energies huj^i, — -\- e^. To allow for a quantitative 
comparison between the strength distributions obtained 
from TDHFB and unperturbed approximations, we com- 
pute the latter using the same time Fourier technique as 
for the TDHFB case. First, we determine the transition 
amplitudes (/iz^l-FjO) and then the time evolution of the 
observable F — (0|-F|0) within the imperturbed approxi- 
mation. The latter reads 



{tiiy\F\0)\^sin{Lo^a)- 



(24) 



Finally, we apply exactly the same procedure as for the 
TDHFB evolution to extract its strength distribution. 
The resulting unperturbed spectrum is represented by a 
dotted line in Fig. [DJa). We see clearly that the effect 
of the residual interaction is to increase the strength on 
the one hand and, on the other hand, to shift down the 
positions of the peaks. 

To get a deeper insight into the nature of the peaks, 
we decompose the response into components of the sin- 
gle particle orbital momentum I using the observables Fi 
defined in Eq. ([^5]) . The response for I — 0,1,2 and 3 
are plotted in Fig. [DJb). These spectra, together with 
the quasiparticle HFB spectrum, allow us to character- 
ize the peaks in terms of dominating two-quasiparticle 
excitations. As one can see in Fig. [2Ib), the first peak 
located at 6.5 MeV is associated to the I = component 
of F. It corresponds mainly to a pair transfer toward the 
almost empty 2si/2 orbitals. 

The next peak, located at 14.3 MeV, is mainly a mix- 
ture of two contributions: the transfer of a pair towards 
1^3/2 orbitals and the removal of the lpi/2 occupied neu- 
trons indicated by a dotted arrow in Fig.[2][a). The fact 
that these two modes have the same energy is fortuitous. 
(This is also the case at the unperturbed level.) As we 
will see later, they are well separated in the other Oxy- 
gen isotopes. We also see in Fig. [2jb) that there is a 
I = 3 contribution to this peak due to a coupling to /2 
orbitals in the continuum. 

Let us now focus on the group of peaks at higher en- 
ergies. As one can see in Fig. [Sfb), they are mostly pop- 
ulated by Z = 1 and 3 components. In fact, the peaks 
between 20 and 24 MeV are mainly associated to the 
excitation of f-j/2 quasiparticle resonant states while the 
peak at 24.7 MeV, indicated by a solid arrow in Fig.[2fa), 
corresponds to the deep hole 1^3/2 state. Except for the 
latter contribution, which is due to the removal of two 
occupied neutrons, these peaks belong to the GPV [s^ ]. 
Indeed, they correspond to excitations of resonant states 
belonging to the next major shell and the enhancement of 
the strength as compared to the unperturbed spectrum 
is a sign of their collectivity pLSj . 
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FIG. 3: Strength distributions of the two-neutron transfer op- 
erator F for i*'20,22q f^_^ arbitrary units with the same scale 
on each plot). TDHFB results (solid lines) and the unper- 
turbed approximation (dotted lines) are shown. The arrows 
indicate the lp3/2 (solid) and lpi/2 (dotted) deep hole states 
and the 1^5/2 (dashed) pair removal. The filled regions cor- 
respond to GPV candidates (see text). 



B. neutron-rich Oxygen isotopes 

In addition to ^^O, we also studied two-neutron pair 
transfer in 20,22q jj^jigj ^jjg spectra are shown in Fig. [3] 
while the energies and most important quasiparticle con- 
tributions to the main peaks are summarized in table [J 
Comparing strength distributions obtained from TDHFB 
(solid lines) with the unperturbed approximation (dotted 
lines) in Fig. [3] leads to the same conclusions for all iso- 
topes, i.e., an increase of the strength and a lowering 
of the peak energies due to the TDHFB residual interac- 
tion. We also see in Fig. [3] that the energies of the IP3/2 
and lj5i/2 deep-hole states increase with the number of 
neutrons. The 1^3/2 peak in ^^O is located outside of 
the figure at 31.3 MeV. The occupied single particle or- 
bitals are indeed deeper as compared to the Fermi level 
when increasing the neutron number, corresponding to 
higher quasiparticle energies. We also note that transi- 



TABLE I: Energies and main quasiparticle contributions of 
the most important peaks appearing in the strength distri- 
bution of the two-neutron pair transfer operator F extracted 
from TDHFB calculations for various Oxygen isotopes. Num- 
bers in parentheses are the centroids of the continuum-QRPA 
energies of Ref . [5^ . Labels in brackets indicate two-neutron 
removal contributions. 
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tions associated to the 1^5/2 orbital appear only in the 
^^O spectrum. This is due to the fact that the operator F 
leaves the strongly paired levels almost unchanged [isj . 
and then no significant strength is associated to addition 
and removal of nucleons into a partially occupied level at 
the Fermi surface. This is not the case in ^^O where the 
1^5/2 single particle orbital is almost fully occupied. 

Last but not least, we see that the GPV, indicated by a 
filled region, is present in the three isotopes with similar 
amplitudes. In all cases, the most important contribu- 
tions to the GPV are the excitation of /7/2 quasiparticle 
resonant states. In the present calculations, continuum 
states are discretized due to the finite size of the box, 
inducing a fragmentation of the GPV. We also note that 
the GPV and the other two-neutron additional modes, 
contrary to the removal ones, have a decreasing energy 
with the neutron number. This is due to the fact that 
the Fermi level is less deep for neutron-rich nuclei, which 
decreases the quasiparticle energy of states with small or 
zero occupation number. 

Let us now compare the energies predicted by the 
present TDHFB calculations with the continuum-QRPA 
results of Khan et al. [52^ , also reported in table [D The 
latter have been computed for two-neutron additional 
modes only. Compared to QRPA results, TDHFB cal- 
culations globally predict slightly higher energies when 
going to more neutron-rich nuclei. In the GPV region, 
the centroid energies are 0.5, 1 and 2 MeV higher with 
TDHFB for i8,20,22q respectively. However, this overall 
agreement can be considered as good in regard to the dif- 
ferences between the two approaches. Both calculations 
use the SLy4 Skyrme functional, but with different pair- 
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ing schemes. In the present work, we use a mixed vohime- 
surface effective couphng. In addition, our calculations 
are performed in wide quasiparticle energy and angular 
momentum windows, with cutoff values Ec — 80 MeV 
and jc — 19/2 respectively, while the QRPA calculations 
have been performed in smaller windows {Ec = 50 MeV 
and jc = 9/2) with a surface type pairing functional. 
The parameters of the latter have been determined using 
a different prescription than ours [l^ , in particular to re- 
produce the trend of the experimental gap in neutron-rich 
Oxygen isotopes. Another possible source of discrepan- 
cies is the fact that the QRPA calculations of Ref. [52] do 
not take into account the Coulomb and spin-orbit parts 
of the residual interaction whereas the TDHFB approach 
uses the same EDF as the underlying HFB field. This 
assumption may induce a slight shift in the energy of 
collective modes [UllKlil. 



C. Calcium isotopes 

In this section we discuss two-neutron pair transfer on 
42,44,46(^g^_ We have plotted in Fig. d] the corresponding 
TDHFB strength distributions (solid lines). The spectra 
are roughly similar for the three isotopes. They exhibit 
several discrete transitions to bound states together with 
excitations of resonant two-quasiparticle states. The en- 
ergy threshold for the latter can be estimated by twice 
the Fermi energy, i.e., 2Ep = 20.6, 19.6 and 17.8 MeV 
for *^Ca, ^'^Ga. and "^^Ca respectively. 

We performed the same analysis as for the Oxygen iso- 
topes, i.e., we decomposed the strength distribution in 
^-components which, together with the HFB quasiparticle 
spectra, helped us assign the main quasiparticle contri- 
butions to each peak. A summary of the results is given 
in tablclTTlwherc the transitions associated to the removal 
of two neutrons arc indicated in brackets. 

For the three isotopes, the lowest mode is interpreted 
in terms of the addition of a neutron pair in the 2p^/2 
orbitals. In "^^Ca, an additional I = 2 quasiparticle com- 
ponent contributes to this mode and corresponds to the 
removal of a 1^3/2 neutron pair. As for Oxygen isotopes, 
the appearance of removal modes (in brackets in table HI]) 
together with additional modes at the same energy is for- 
tuitous. In the Calcium isotopes, the removal modes are 
built of neutrons from the s — d shell, the major shell be- 
low the Fermi energy. As expected, one finally notes that 
energies of the removal (resp. additional) modes increase 
(decrease) with the neutron number. Although 2 179/2 
quasiparticles excitations are forbidden below the 2Ef 
threshold in the unperturbed spectrum, it gets mixed to 
the last bound 2 (75/2) quasiparticles excitation because 
of the residual interaction. 

We have also plotted the unperturbed spectra in Fig. [3] 
(dotted lines). As in the case of Oxygen isotopes, we ob- 
serve that the TDHFB residual interaction lowers the en- 
ergies on the one hand and increases the strength on the 
other hand, though this second effect is less pronounced 
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FIG. 4: Strength distributions of the two-neutron pair trans- 
fer operator F for ''^■^^''''^Ca (in arbitrary units witii the same 
scale on each plot). TDHFB results (solid lines) and the 
unperturbed approximation (dotted lines) are shown. The 
arrows indicate the ld5/2 (solid), 2si/2 (dashed) and 1^3/2 
(dotted) two deep-hole states. 



in the high energy part of the spectra. Here, the pairing 
residual interaction is not strong enough to gather high 
energy states together and to generate a well identified 
collective pairing vibration in the continuum. In a pure 
independent particles shell model picture, the last occu- 
pied neutron level of these Calcium isotopes is the I/7/2 
orbital. Then, one expects the GPV to be built mainly 
on low-lying gg/2 resonant quasiparticle states. However, 
as we clearly see in Fig. 21 the strength associated to 
the two-quasiparticle excitation of 59/2 levels, located at 
23.4, 22.1 and 20.8 MeV for ^SQa, ^-^Ca and "^^Ca respec- 
tively, is only slightly enhanced by the TDHFB residual 
interaction. We checked that employing other parameter- 
izations of the spatial distribution /(r) of the excitation 
operator in Eq. ([2T|) does not alter these conclusions. We 
note that this lack of collectivity of the gg/2 has already 
been observed in hole pairing giant resonances studies 
within a more schematic formalism [77i | . 
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TABLE II: Energies and main quasiparticle contributions of 
the most important peaks appearing in the strength distri- 
bution of the two-neutron pair transfer operator F extracted 
from TDHFB calculations for Calcium isotopes. Labels in 
brackets indicate two-neutron removal contributions. 
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these differences. 

In addition, there is room for a better treatment of the 
continuum, for instance in the spirit of continuum-QRPA 
calculations, but such an improvement is not straight- 
forward. Indeed, one cannot extrapolate the absorbing 
boundary conditions used in TDHF calculations to the 
TDHFB case because of the delocalized upper compo- 
nents of the Bogoliubov spinors. 

Finally, TDHFB calculations are much more demand- 
ing in terms of computational time than standard TDHF 
calculations, by about two orders of magnitude in the one 
dimensional case. However, thanks to recent increase of 
computational power, some of the standard TDHF ap- 
plications to nuclear structure and reactions should be 
repeated with the inclusion of dynamical pairing correla- 
tions in the framework of TDHFB. 
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VI. CONCLUSION 

We solved the TDHFB equation in coordinate space 
with spherical symmetry for the evolution of a single nu- 
cleus in an external field. For the normal part of the en- 
ergy density functional, we used the SLy4 Skyrme func- 
tional. For its pairing part, we chose a local density de- 
pendent functional. Special care has been taken regard- 
ing the convergence of the static HFB solutions and the 
energy and particle number conservations in the TDHFB 
calculations. 

As a first application, we studied 0"*" pairing modes ex- 
cited by a two-neutron pair transfer type operator. The 
linear response theory has been used to compute the 
strength distributions of this operator in i8,20,22q 
42,44,46(-ig^ nuclei. Both transitions to bound states and 
to the continuum are observed in all nuclei. In particular, 
the GPV is observed in all Oxygen isotopes whereas no 
significant enhancement of the strength due to dynami- 
cal pairing correlations appears in the continuum of the 
studied Calcium isotopes. In the latter, the ggy2 quasi- 
particle excitations are not collective enough to generate 
a GPV. 

A detailed comparison with previous QRPA calcu- 
lations have been performed in the Oxygen isotopes. 
Though there is a good agreement for the most stable 
isotope, we find slightly higher energies for the pairing 
vibrations when going to more neutron-rich nuclei. Dif- 
ferent pairing schemes and implementations of the resid- 
ual interaction in both calculations are invoked to explain 



APPENDIX A: DENSITIES AND FIELDS IN 
SPHERICAL SYMMETRY 

Each energy term entering Eq. ([5]) can be written as a 
spatial integral of an energy density, i.e., Z = /dr H(r), 
which depends only on local densities and currents. For 
a spherically symmetric system, the densities entering 
the SLy4 and local density dependent pairing functionals 
are the radial part of the local particle Pq{j\ anomalous 
Pq(r), kinetic Tg(r), matter current \q{y) and spin-orbit 
current Jg(r) (both oriented along the radial unit vector 
e^) densities of isospin q. Introducing the notation a = 
{n, /, J, q), these densities write 

p,(r) - ^if,(r)|i;„(r)p, 

nlj 

Pq{^) = -^Kj{r)vl(r)ua{r), 

nlj 



r,(r) - Y.KAr) 

nlj 



I- - J ) 



nlj 



nlj 



\va{r)\'^er 
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where Kj{r) = '^JJ^l . We also define isoscalar densities 
as the sum of proton and neutron densities, e.g., /Oo(r) 
Pp{r) + jO„(r). Omitting the dependence in r to simphfy 
the notations, the different parts of the functional can be 
written 



'Hsk = 

Ti-Coul = 



k—n,p,0 



~Cl{pkTk~3l)+C'pkV.3k } 



rrdir 32/3 



q=p,n 
1/3 



/"IT 

^n,p 



Xl 



xi + - 



X2 



16 
16'^^ 



r 

2 
1 

' 2 



The fields entering Eq. write 

d d d 

hi,,{r) = -_Af,(r)-+/,(r)- + l^,,(r) 
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Po(?-) 



with 



where V/"" (r) — ^ J dr'j^^ is the direct Coulomb 

field. The factor (1-1 /A) in the kinetic part is the so- 
called one-body center of mass correction. The term 
ensures Galilean invariance [TSj . 

The coefficients C are related to the usual Skyrme co- 
efficients by (see, e.g., [ssj ) 



^0 



^n,p 



to 
2 



(i + f), 

to f ^1 

'Yr" + 2 



^fi + ^ 

12 V 2 



^3 

12 



2:3 



4 



X2 

2 



2m 



+ Qpo + C^jPg 



I, = 2*(Qio + C;j,).e, 



dr 



U„ = 



r 



E [2 {c^k+cCpojPk 

k=q,0 

+aCCpr\pl + pl) + ClTk 



7-1 

5 Po l\~ \2 

-i7^(|Pp| 



[1] J. R. Stone and P.-G. Reinhard, Prog. Part. Nucl. Phys. 

58, 587 (2007). 
[2] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. 

Phys. 75, 121 (2003). 
[3] D. Vautherin and D. M. Brink, Phys. Rev. C 5, 626 

(1972). 

[4] M. Bender, G. F. Bertsch, and P.-H. Heenen, Phys. Rev. 

Lett. 94, 102503 (2005). 
[5] M. Bender, G. F. Bertsch, and P.-H. Heenen, Phys. Rev. 

C 73, 034322 (2006). 
[6] A. S. Umar and V. E. Oberacker, Phys. Rev. C 73, 

054607 (2006). 

[7] C. Simenel and B. Avez, Int. J. Mod. Phys. E 17, 31 
(2008). 

[8] D. R. Hartree, Proc. Cambridge Philos. Soc. 24, 89 
(1928). 

[9] V. Fock, Z. Phys. A 61, 126 (1930). 
[10] N. Bogoliubov, Sov. Phys. JETP 7, 41 (1958). 
[11] J. Engel, M. Bender, J. Dobaczewski, W. Nazarewicz, 

and R. Surman, Phys. Rev. C 60, 014302 (1999). 
[12] M. Bender, J. Dobaczewski, J. Engeb and 



W. Nazarewicz, Phys. Rev. C 65, 054322 (2002). 
[13] E. Khan, N. Sandulescu, M. Grasso, and N. Van Giai, 

Phys. Rev. C 66, 024309 (2002). 
[14] S. Fracasso and G. Colo, Phys. Rev. C 72, 064310 (2005). 
[15] J. Terasaki, J. Engel, M. Bender, J. Dobaczewski, 

W. Nazarewicz, and M. Stoitsov, Phys. Rev. C 71, 

034310 (2005). 

[16] J. Terasaki and J. Engel, Phys. Rev. C 74, 044301 (2006). 
[17] S. Peru and H. Goutte, Phys. Rev. C 77, 044313 (2008). 
[18] D. Lacroix, S. Ayik, and P. Chomaz, Prog. Part. Nucl. 

Phys. 52, 497 (2004). 
[19] T. Nakatsukasa and K. Yabana, Phys. Rev. C 71, 024301 

(2005). 

[20] G. Colo and N. Van Giai, Nucl. Phys. A731, 15 (2004). 
[21] S. Peru, J.-F. Berger, and P.-F. Bortignon, Eur. Phys. J. 

A 26, 25 (2005). 
[22] T. Sil, S. Shlomo, B. K. Agrawal, and P. G. Reinhard, 

Phys. Rev. C 73, 034316 (2006). 
[23] M. Holland, J. Park, and R. Walser, Phys. Rev. Lett. 86, 

1915 (2001). 

[24] H. Buljan, M. Segev, and A. Vardi, Phys. Rev. Lett. 95, 



11 



180401 (2005). [52 
J. W. Ncgclo, S. E. Koonin, P. MoUcr, J. R. Nix, and 
A. J. Sicrk, Phys. Rev. C 17, 1098 (1978). [53 
R. Y. Cusson and H. W. Meldner, Phys. Rev. Lett. 42, 
694 (1979). [54 
R. Y. Cusson, J. A. Maruhn, and H. Stocker, Zeit. Phys. 
A 294, 257 (1980). [55 
J. Blocki and H. Flocard, Nucl. Phys. A273, 45 (1976). [56 
J. Dcchargo and D. Gogny, Phys. Rev. C 21, 1568 (1980). 
Y. Hashimoto and K. Nodeki, arXiv : 0707.3083 [nucl-th] [57 
(2007); talk given at the FIDIPRO-JSPS Workshop on 
Energy Density Punctionals in Nuclei, 25-27 Oct. 2007, [58 
Keurusselka, Jyvaskyla (Finland). 
P. A. M. Dirac, Proc. Cambridge Philos. Soc. 26, 376 [59 
(1930). 

T. Skyrme, Phil. Mag. 1, 1043 (1956). [60 
P. Bonche, S. Koonin, and J. W. Negele, Phys. Rev. C 
13, 1226 (1976). 

J. W. Negele, Rev. Mod. Phys. 54, 913 (1982). [61 
C. Simcnol, P. Chomaz, and G. de Prance, Phys. Rev. 
Lett. 86, 2971 (2001). [62 
C. Simenel, P. Chomaz, and G. de Prance, Phys. Rev. [63 
Lett. 93, 102701 (2004). 
J. A. Maxuhn, P.-G. Reinhaxd, P. D. Stevenson, and [64] 
M. R. Strayer, Phys. Rev. C 74, 027601 (2006). [65 
C. Simenel, P. Chomaz, and G. de Prance, Phys. Rev. C 
76, 024609 (2007). [66 

C. Simenel and P. Chomaz, Phys. Rev. C 68, 024302 
(2003). [67 
A. S. Umar and V. E. Oberacker, Phys. Rev. C 71, 
034314 (2005). [68 
,J. A. Maruhn, P.-G. Reinhard, P. D. Stevenson, J. R. [69 
Stone, and M. R. Strayer, Phys. Rev. C 71, 064328 

(2005) . 

P.-G. Reinhard, L. Guo, and J. A. Maruhn, Eur. Phys. [70 
J. A 32, 19 (2007). [71 
A. Bohr and B. Mottelson, Nuclear Structure (2 vol., [72 
W.A. Benjamin, Inc., 1975). 

P. Ring and P. Schuck, The Nuclear Many-Body Problem [73 

(Springer Verlag, 1980). 

D. Res and R. Brogha, Nucl. Phys. A80, 289 (1966). [74 
C. H. Dasso and M. I. Gallardo, Phys. Rev. C 74, 014307 

(2006) . [75 
R. Broglia and D. Bcs, Phys. Lett. B69, 129 (1977). 
G. Ripka and R. Padjen, Nucl. Phys. A132, 489 (1969). [76 
R. Broglia, O. Hansen, and C. Riedel, Advances in Nu- 
clear Physics 6, 287 (1973). [77 
L. Fortunato, W. von Oertzen, H. M. Sofia, and A. Vit- 

turi, Eur. Phys. Jour. A 14, 37 (2002). [78 
W. von Oertzen and A. Vitturi, Rep. Prog. Phys. 64, 
1247 (2001). 



E. Khan, N. Sandulescu, N. Van Giai, and M. Grasso, 
Phys. Rev. C 69, 014314 (2004). 

J.-P. Blaizot and G. Ripka, Quantum Theory of Finite 
Systems (MIT Press, 1986). 

P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 
(1964). 

W. Kohn and L. Sham, Phys. Rev. 140, A1133 (1965). 
E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 

(1984). 

E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and 
R. Schaeffer, Nucl. Phys. A635, 231 (1998). 
P. Bonche, H. Flocard, and P.-H. Heenen, Nucl. Phys. 
A467, 115 (1987). 

J. W. Negele and D. Vautherin, Nucl. Phys. A207, 298 
(1973). 

J. Dobaczewski, W. Nazarewicz, T. R. Werner, J. F. 
Berger, C. R. Chinn, and J. Decharge, Phys. Rev. C 53, 
2809 (1996). 

G. Bruun, Y. Castin, R. Dum, and K. Burnett, Eur. 
Phys. J. D 7, 433 (1999). 

A. Bulgac and Y. Yu, Phys. Rev. Lett. 88, 042504 (2002). 

J. Dobaczewski, W. Nazarewicz, and P.-G. Reinhard, 

Nucl. Phys. A693, 361 (2001). 

A. Bulgac, Phys. Rev. C 41, 2333 (1990). 

J. Dobaczewski, H. Flocard, and J. Treiner, Nucl. Phys. 

A422, 103 (1984). 

K. Bennaceur and J. Dobaczewski, Comp. Phys. Com. 
168, 96 (2005). 

H. Flocard, S. E. Koonin, and M. S. Weiss, Phys. Rev. 
C 17, 1682 (1978). 

D. Baye, J. Phys. B 28, 4399 (1995). 

P.-G. Reinhard, P. D. Stevenson, D. Almehed, J. A. 

Maruhn, and M. R. Strayer, Phys. Rev. E 73, 036709 

(2006). 

J. Blocki and H. Flocard, Phys. Lett. B85, 163 (1979). 
S. Stringari and D. Vautherin, Phys. Lett. B88, 1 (1979). 
A. S. Umar and M. R. Strayer, Phys. Lett. B171, 353 

(1986) . 

P. Chomaz and N. Van Giai, Phys. Lett. B189, 375 

(1987) . 

M. V. C. R. Chinn, A. S. Umar, and M. R. Strayer, Phys. 
Rep. 264, 107 (1996). 

D. Almehed and P. D. Stevenson, J. Phys. G 31, S1819 

(2005). 

P. Stevenson, D. Almehed, P.-G. Reinhard, and 
J. Maruhn, Nucl. Phys. A788, 343 (2007). 
M. Herzog, R. Liotta, and T. Verste, Phys. Lett. B165, 
35 (1985). 

Y. Engel, D. Brink, K. Goeke, S. Krieger, and D. Vau- 
therin, Nucl. Phys. A249, 215 (1975). 



